import numpy as N
import math as M
from Interior_lightcone import get_solution_region_I
#from Vrai_vers_faux import get_solution_region_II
from Collage import Collage
import spline
import matplotlib.pyplot as plt
import scipy
import scipy.optimize

epsilon = 0.01 # Meme valeur de Match_b_col

# Region II

#tvec_II,yvec_II,f = get_solution_region_II()
tvec_II,D_phi_II,phi_II,D_b,b,f,phi0 = Collage(epsilon)
#b   = yvec_II[:,3]
#D_b = yvec_II[:,2]
#phi_II   = yvec_II[:,1]
#D_phi_II = yvec_II[:,0]
b_func   = spline.Spline(tvec_II,b)
D_b_func = spline.Spline(tvec_II,D_b)
phi_II_func   = spline.Spline(tvec_II,phi_II)
D_phi_II_func = spline.Spline(tvec_II,D_phi_II)
f_func = spline.Spline(tvec_II,f)

# Region I

t_start = 0.005 # Start time of numerical integration
t_final = 3000. # "Duration" of the proccess
pas = 0.005 # Stepsize

tvec_I,yvec_I,g = get_solution_region_I(t_start,t_final,phi0,pas,epsilon)
a   = yvec_I[:,3] # I call it a as well, but I have to change it to avoid misunderstandings
D_a = yvec_I[:,2] # I call it a as well, but I have to change it to avoid misunderstandings
phi_I   = yvec_I[:,1]
D_phi_I = yvec_I[:,0]
a_func   = spline.Spline(tvec_I,a) # I call it a as well, but I have to change it to avoid misunderstandings
D_a_func = spline.Spline(tvec_I,D_a) # I call it a as well, but I have to change it to avoid misunderstandings
phi_I_func   = spline.Spline(tvec_I,phi_I)
D_phi_I_func = spline.Spline(tvec_I,D_phi_I)
g_func = spline.Spline(tvec_I,g)

#def tau(u,v): # Expression (54)
#  aux = M.sqrt(u*v)*M.exp((f_func(u)+f_func(v))/2.) # Premiere estimation
#  result = aux
#  for i in range(100): # Prend beaucoup de temps, et sa convergence est de plus en plus lente pour u,v grands
#	result = aux*M.exp(-g_func(result)) # C'est l'iteration afin d'ameliorer l'estimation
#  return (result)

def tau(u,v): # Expression (54)
  aux = M.sqrt(u*v)*M.exp((f_func(u)+f_func(v))/2.) # Premiere estimation
  def func(tau):
    return (tau*M.exp(g_func(tau)) - M.sqrt(u*v)*M.exp((f_func(u)+f_func(v))/2.))
  return (scipy.optimize.newton(func, aux))

#def sig(u,v): # Expression (54) pour sigma
#  aux = M.sqrt(u*v)*M.exp((f_func(u)+f_func(v))/2.) # Premiere estimation
#  result = aux
#  for i in range(30): # Assez bonne convergence
#  	result = aux*M.exp(-f_func(result))
#  return (result)

def sig(u,v): # Expression (54) pour sigma
  aux = M.sqrt(u*v) # Premiere estimation
  def func(sig):
    return (sig*M.exp(f_func(sig)) - M.sqrt(u*v)*M.exp((f_func(u)+f_func(v))/2.))
  return (scipy.optimize.newton(func, aux, maxiter=100))

def a_exact(u,v): # Racine de l'expression (62) - cette expression verifie f^2(u,-u)=1.
  if (u<0.): # Reflexion
	aux = u
	u = v
	v = aux
  if (v>0.):
  	result=M.log(a_func(tau(u,v))**2/(b_func(u)*b_func(v)))/2.
  else:
	v=-v # On est sur la region II, a l'exterieur du cone de lumiere, au cote droit (u>0,v<0 avec -v<u)
	result=M.log(b_func(sig(u,v))**2/(b_func(u)*b_func(v)))/2.
  return(result)

def a_u_exact(u,v): # J'ai fait les derivations ; afin de vous aider, sachez que dtau(u,v)/du=a(tau(u,v))/(2b(u)) et dtau(u,v)/dv=a(tau(u,v))/(2b(v)) (calcul base sur (54) et (60))
  if (u<0.): # Reflexion
	aux = u
	u = v
	v = aux
	return(a_v_exact(u,v))
  if (v>0.):
  	result= (D_a_func(tau(u,v)) - D_b_func(u))/(2.*b_func(u))
  else:
	v=-v
	result= (D_b_func(sig(u,v)) - D_b_func(u))/(2.*b_func(u))
  return(result)

def a_v_exact(u,v):
  if (u<0.): # Reflexion
	aux = u
	u = v
	v = aux
	return(a_u_exact(u,v))
  if (v>0.):
  	result= (D_a_func(tau(u,v)) - D_b_func(v))/(2.*b_func(v))
  else:
	v=-v
	result=-(D_b_func(sig(u,v)) - D_b_func(v))/(2.*b_func(v))
  return(result)

def r_exact(u,v): # Racine de l'expression (63) ; on utilise (47) et (48) pour avoir les expressions de t et s - cette expression verifie r(u,-u)=0.
  if (u<0.): # Reflexion
	aux = u
	u = v
	v = aux
  if (v>0.):
  	result= a_func(tau(u,v))*M.cosh((M.log(u/v)+f_func(u)-f_func(v))/2.)
  else:
	v=-v
  	result= b_func(sig(u,v))*M.sinh((M.log(u/v)+f_func(u)-f_func(v))/2.)
  return(result)

def r_u_exact(u,v): # J'ai fait les derivations ; afin de vous aider, sachez que dt+(tau+)/du=ds+(u)/du=1/b(u) et dt-(tau-)/dv=ds-(v)/dv=1/b(v)
  if (u<0.): # Reflexion
	aux = u
	u = v
	v = aux
	return(r_v_exact(u,v))
  if (v>0.):
  	result= a_func(tau(u,v))*(D_a_func(tau(u,v))*M.cosh((M.log(u/v)+f_func(u)-f_func(v))/2.)+M.sinh((M.log(u/v)+f_func(u)-f_func(v))/2.))/(2.*b_func(u))
  else:
	v=-v
	result= b_func(sig(u,v))*(D_b_func(sig(u,v))*M.sinh((M.log(u/v)+f_func(u)-f_func(v))/2.)+M.cosh((M.log(u/v)+f_func(u)-f_func(v))/2.))/(2.*b_func(u))
  return(result)

def r_v_exact(u,v):
  if (u<0.): # Reflexion
	aux = u
	u = v
	v = aux
	return(r_u_exact(u,v))
  if (v>0.):
  	result= a_func(tau(u,v))*(D_a_func(tau(u,v))*M.cosh((M.log(u/v)+f_func(u)-f_func(v))/2.)-M.sinh((M.log(u/v)+f_func(u)-f_func(v))/2.))/(2.*b_func(v))
  else:
	v=-v
	result=-b_func(sig(u,v))*(D_b_func(sig(u,v))*M.sinh((M.log(u/v)+f_func(u)-f_func(v))/2.)-M.cosh((M.log(u/v)+f_func(u)-f_func(v))/2.))/(2.*b_func(v))
  return(result)

def phi_exact(u,v):
  if (u<0.): # Reflexion
	aux = u
	u = v
	v = aux
  if (v>0.):
  	result=phi_I_func(tau(u,v))
  else:
	v=-v
	result=phi_II_func(sig(u,v))
  return(result)

def phi_u_exact(u,v):
  if (u<0.): # Reflexion
	aux = u
	u = v
	v = aux
	return(phi_v_exact(u,v))
  if (v>0.):
  	result= D_phi_I_func(tau(u,v))*a_func(tau(u,v))/(2.*b_func(u))
  else:
	v=-v
	result= D_phi_II_func(sig(u,v))*b_func(sig(u,v))/(2.*b_func(u))
  return(result)

def phi_v_exact(u,v):
  if (u<0.): # Reflexion
	aux = u
	u = v
	v = aux
	return(phi_u_exact(u,v))
  if (v>0.):
  	result= D_phi_I_func(tau(u,v))*a_func(tau(u,v))/(2.*b_func(v))
  else:
	v=-v
	result=-D_phi_II_func(sig(u,v))*b_func(sig(u,v))/(2.*b_func(v))
  return(result)
